function [A, A0] = sllinreg(X, Y, varargin)
%SLLINREG Performs linear regression on data
%
% [ Syntax ]
%   - [a, a0] = sllinreg(X, y, ...)
%   - [A, A0] = sllinreg(X, Y, ...)
%
% [ Argument ]
%   - X:        The design matrix of regression (m x n)
%   - y:        The response vector (m x 1)
%   - Y:        The response matrix (m x k)
%   
%   - a:        The solved coefficient vector (n x 1)
%   - a0:       The constant offset (1 x 1)
%   - A:        The coeffient matrix (n x k)
%   - A0:       The vector of constant offsets (1 x k)
%
% [ Description ]
%   - [a0, a] = sllinreg(X, y, ...) solves the linear regression problem, 
%     which can be formalized as follows
%       $ minimize (1/2) ||y - Xa - a0||^2 + lambda^T R(a) $
%
%     Let X be an m x n matrix, and y be an m x 1 vector. Then it solves
%     the linear regression problem and results in a coefficient vector
%     a of size n x 1 and a constant scalar a0.
%
%     By default the problem is unregularized. However, you can specify 
%     how to regularize the coefficients by options.
%
%     Here are available options
%     \{:
%           - lambda:       The the regularization coefficient, which can 
%                           be a scalar which imposes uniform
%                           regularization on all entries of a, or a vector
%                           of size n x 1 to impose different
%                           regularization levels to different entries.
%                           all values in lambda should be positive.
%                           (By default, lambda = 0, which implies that
%                           no regularization is imposed)
%                           
%                           Note that lambda can be a vector only when reg
%                           is L1 or L2.
%
%           - reg:          The form of regularization. (default = 'L2')
%                           \{:
%                               - L2:   The squared L2 norm regularization
%                                       $ R(a) = (1/2) ||a||_2^2 $
%                               - L1:   The L1 norm regularization
%                                       $ R(a) = |a|_1 $
%                               - GTik: Generic Tikhonov regularization
%                                       $ R(a) = (1/2) a^T Q ^a $                                                                    
%                           \:}                         
%
%           - use_offset:   Whether to enable offset constant. The default
%                           is true. If you set it to be false, then a0 is
%                           kept at zero.
%
%           - W:            The weights of samples. By default, W = 1.
%                           Depending on its form, the first term of the 
%                           objective function will change.
%                           If W is a scalar, then it is W * ||err||^2;
%                           If W is a vector of size m x 1, then it is
%                           \sum w_i err_i
%                           If W is a squared matrix, then it is 
%                           err^T * W * err. 
%                           Please note that when W is scalar or vector,
%                           all its values should be positive, when it is
%                           a matrix, it should be positive definite.
%
%           - Q:            The matrix used in Generic Tikhonov
%                           regularization.By default, Q = []. If you
%                           set the option reg to 'GTik', you need to
%                           specify a non-empty n x n positive
%                           semi-definite matrix for Q.
%
%           - solver:       If the problem is unregularized or regularized
%                           with L2 or Generalized Tikhonov, it involves 
%                           solving a linear equation (A' A)X = A' y. 
%                           By default, we use standard way: (A'A) \ (A'y).
%                           However, based on the positive definiteness of
%                           A'A, some different methods may be more
%                           efficient. 
%                           \{:
%                               - []:       The standard way using mldivide
%                               - pcg:      Using Preconditioned Conjugate
%                                           Gradient method
%                               - minres:   Using Minimal Residual Method
%                               - qmr:      Using Quasi-Minimal Residual
%                                           Method
%                               - symmlq:   Using Symmetric LQ Method
%                           \:}
%                           You can also spesify your solver in form of a
%                           function handle with syntax X = f(A, B) to
%                           solve equation system AX = B, where A would be
%                           positive semidefinite. 
%
%                           If you are using L1 regularizer, the problem
%                           will be casted to a linearly constraint
%                           quadratic programming problem, and the function
%                           will invoke a QP solver.  
%                           By default, it calls the quadprog. You can also
%                           specify your own QP solver by setting this
%                           option as a function handle with the syntax 
%                           x = f(H, f, A, b).
%                           
%     \:}
%
%   - [A0, A] = sllinreg(X, Y, ...) solves k linear regression problem
%     sharing the same design matrix in batch. 
%
% [ History ]
%   - Created by Dahua Lin, on Oct 31, 2007
%

%% parse and verify input arguments

error(nargchk(2, inf, nargin));

% X and Y
assert(isnumeric(X) && ndims(X) == 2, ...
    'sltoolbox:sllinreg:invalidarg', ...
    'X should be a 2D numeric matrix.');
[m, n] = size(X);

assert(isnumeric(Y) && ndims(Y) == 2, ...
    'sltoolbox:sllinreg:invalidarg', ...
    'y or Y should be a numeric vector or matrix.');

assert(size(Y, 1) == m, ...
    'sltoolbox:sllinreg:invalidarg', ...
    'The first dimension of y or Y should be m.');
k = size(Y, 2);

% options
opts = struct( ...
    'lambda', 0, ...
    'reg', 'L2', ...
    'use_offset', true, ...
    'W', 1, ...
    'Q', [], ...
    'x0', [], ...
    'solver', []);
opts = setopts(opts, varargin{:});

lambda = opts.lambda;
assert(isnumeric(lambda) && (isscalar(lambda) || isequal(size(lambda), [n 1])), ...
    'sltoolbox:sllinreg:invalidopt', ...
    'lambda should be either a scalar or a n x 1 vector.');

if all(lambda == 0)
    do_reg = false;
else
    do_reg = true;
    assert(all(lambda >= 0), ...
        'sltoolbox:sllinreg:invalidopt', ...
        'all values in lambda should be nonnegative.');
end

if do_reg
    reg = opts.reg;
    assert(ischar(reg), ...
        'sltoolbox:sllinreg:invalidopt', ...
        'The option reg should be a char string.');     
    
    switch reg
        case {'L2', 'GTik'}
            is_qp = false;
        case 'L1'
            is_qp = true;
        otherwise
            error('sltoolbox:sllinreg:invalidopt', ...
                'Unknown regularizer form %s', reg);
    end
end

uof = opts.use_offset;
assert(isscalar(uof) && islogical(uof), ...
    'sltoolbox:sllinreg:invalidopt', ...
    'The option use_offset should be a logical scalar.');

W = opts.W;
assert(isnumeric(W), ...
    'sltoolbox:sllinreg:invalidopt', ...
    'The option W should be of numeric type');

if isscalar(W)
    assert(W > 0, 'sltoolbox:sllinreg:invalidopt', ...
        'W should be nonnegative.');
    
elseif isequal(size(W), [m 1])
    assert(all(W) > 0, 'sltoolbox:sllinreg:invalidopt', ...
        'All values in W should be nonnegative.');
    if all(W == W(1))
        W = W(1);
    end
    
elseif isequal(size(W), [m m])
    if isequal(W, W(1) * eye(m, m))
        W = W(1);
    elseif isequal(W, diag(diag(W)))
        W = diag(W);
    end
    
else
    error('sltoolbox:sllinreg:invalidopt', ...
        'The option W should be either a scalar, an m x 1 vector, or an m x m matrix.');
end

if strcmp(opts.reg, 'GTik')
    assert(isscalar(lambda), ...
        'sltoolbox:sllinreg:invalidopt', ...
        'lambda should be a scalar when using GTik');

    Q = opts.Q;
    assert(isnumeric(Q) && isequal(size(Q), [n n]), ...
        'sltoolbox:sllinreg:invalidopt', ...
        'Q should be specified as an n x n matrix for GTik.');
end


if ~is_qp
    if isempty(opts.solver)
        fh_solver = @mldivide;
    elseif ischar(opts.solver)
        switch opts.solver
            case 'std'
                fh_solver = @mldivide;
            case 'pcg'
                fh_solver = @pcg;
            case 'minres'
                fh_solver = @minres;
            case 'qmr'
                fh_solver = @qmr;
            case 'symmlq'
                fh_solver = @symmlq;
            otherwise
                error('sltoolbox:sllinreg:invalidopt', ...
                    'Unknown solver name %s', solver);
        end
        
    elseif isa(opts.solver, 'function_handle')
        fh_solver = opts.solver;
        
    else
        error('sltoolbox:sllinreg:invalidopt', ...
            'solver is invalidly specified');
    end
else
    if isempty(opts.solver)
        solopt = optimset('quadprog');
        solopt = optimset(solopt, 'Display', 'iter', 'LargeScale', 'off');        
        fh_solver = @(H, f, A, b) quadprog(H, f, A, b, [], [], [], [], [], solopt);
        
    elseif isa(opts.solver, 'function_handle');
        fh_solver = opts.solver;
        
    else
        error('sltoolbox:sllinreg:invalidopt', ...
            'solver is invalidly specified');
    end
end



%% Main part

% Unregularized problem

if ~do_reg
    [H, F] = makeHF(X, Y, uof, W);
    
    A = fh_solver(H, F);
    
else
    
    % Regularized problem
        
    switch reg
        
        % L2 
        
        case 'L2'
            [H, F] = makeHF(X, Y, uof, W);
            
            if uof
                diag_inds = 1+(n+2)*(0:(n-1))';            
            else
                diag_inds = 1+(n+1)*(0:(n-1))';
            end
            H(diag_inds) = H(diag_inds) + lambda;
            
            A = fh_solver(H, F);

        % GTik
        
        case 'GTik'
            if lambda ~= 1
                Q = Q * lambda;
            end
                        
            [H, F] = makeHF(X, Y, uof, W);
            
            if uof
                H(1:n, 1:n) = H(1:n, 1:n) + Q;
            else
                H = H + Q;
            end
            
            A = fh_solver(H, F);
            
        % L1
        
        case 'L1'
            [H, F] = makeHF(X, Y, uof, W);
            
            if issparse(H)
                H = blkdiag(H, sparse(n, n));
            else
                H = blkdiag(H, zeros(n, n));
            end            
            
            if isscalar(lambda)
                fu = lambda * ones(n, 1);
            else
                fu = lambda;
            end
            
            if uof
                ds = n + 1;
            else
                ds = n;
            end
            
            A_ie = [-speye(n, ds), -speye(n, n); ...
                     speye(n, ds), -speye(n, n)];
            b_ie = zeros(2 * n, 1);
            
            if k == 1
                f = [-F; fu];
                aa = fh_solver(H, f, A_ie, b_ie);   
                A = aa(1:ds);
            else                                            
                A = zeros(ds, k);

                for i = 1 : k
                    fx = F(:, i);
                    f = [-fx; fu];

                    aa = fh_solver(H, f, A_ie, b_ie);
                    A(:, i) = aa(1:ds);
                end
            end
            
    end % end switch
     
end


% Make the output
   
if uof
    A0 = A(n+1, :);
    A = A(1:n, :);
else
    A0 = zeros(1, size(A, 2));
end


%% Core computation routine

function [H, F] = makeHF(X, Y, uof, W)

m = size(X, 1);

if uof
    Z = [X, ones(m, 1)];
else
    Z = X;
end

if isscalar(W)
    H = Z' * Z;
    F = Z' * Y;
    if W ~= 1
        H = W * H;
        F = W * F;
    end
    
elseif isvector(W)
    H = Z' * bsxfun(@times, W, Z);
    if size(Y, 2) == 1
        F = Z' * (W .* Y);
    else
        F = Z' * bsxfun(@times, W, Y);
    end
        
else    
    H = Z' * W * Z;
    F = Z' * W * Y;
    
end


